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STUDY OF SOME NUMERICAL ARTIFACTS INTERVENING 
IN THE FINITE DIFFERENCES SIMULATIONS OF KdV 
SOLITONS PROPAGATION 


Ion APOSTOL!, Dan Alexandru IORDACHE’, Pier Paolo DELSANTO*, 
Viorica IORDACHE* 


Abstract. As it is well known, the complex simulations are frequently affected by 
numerical artifacts, which reduce drastically their accuracy. Taking into account that the 
numerical artifacts due to the non-linearity of the considered physical problem are less 
studied that those due to the symmetry breaking or to the high logical depth of the used 
schemes, the finite differences (FD) simulations of some (Korteweg-de Vries) solitonic 
pulses were studied. In this aim, both the possible FD discretizations and the artifacts 
corresponding to different values of the simulation uniqueness parameters were studied. 
The possibilities to explain quantitatively the observed numerical artifacts, transforming 
them in numerical phenomena, were also studied. 


Keywords: Complex numerical simulations, numerical artifacts, Korteweg-de Vries solitary 
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1. Introduction 


As it is well-known, we live in a computerized world, our civilization being a 
“civilization of computers”. Since all complex installations and devices are now 
controlled by computers, their corresponding complex simulations can lead to 
some numerical artifacts and — if they are used to control such systems — to major 
failures. Particularly, the Patriot missile failure to stop a Scud missile during the 
Gulf War in 1991 (with disastrous results) [1] and the self-destruction of the 
European Space Agency’s Ariane 5 rocket, at 37 seconds after its launch [2], were 
both assigned to some computer errors [3] and their associated numerical artifacts. 
If the mechanisms of the computer artifacts installation can be quantitatively 
described, these artifacts are called « numerical phenomena » [4]. 


A typical source of numerical artifacts corresponds to the Finite Difference (FD) 
simulations of the solitary waves propagation, due to the nonlinear character of 
their equations. Taking into account the major applications of the more than 100 
types of different solitary waves [5] in all Physics fields: Cosmology and 
Astrophysics, Mechanics and Fluid Mechanics, Molecular chains, Condensed 
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matter Physics (involving the applications for the information transmission 
through optical fibers [6], the Josephson (superconductive) devices, etc.), Plasma 
Physics, Nuclear Physics, applications in Medicine [7] etc. 


A detailed analysis of numerical artifacts intervening in the FD simulations of 
pulses propagation [4], [8], pointed out that — for the complex simulations due to: 
a) the symmetry breaking and to: b) the high logical depth of the involved 
computer programs - their basic types.are.those of: (i) instability, (1i) divergence 
or pseudo-convergence (even more dangerous, because it produces some 
misleading results), (ii1) distortions, etc. 


The main goal of this work is to point out the specific numerical artifacts 
corresponding to the FD simulations associated to pulses propagation described 
by nonlinear equations. 


2. Studied FD discretizations of the KdV equation 


Starting from the one-dimensional (1D) propagation equation [9] of the 
Korteweg-de Vries (KdV) solitons, in terms of their non-linearity n and dispersion 
d coefficients, respectively, and of the unperturbed speed (in absence of the 
nonlinearity and dispersion properties of the propagation medium) V,,: 

Y=—Voo: yy y =a (1) 


the studied discretization schemes are built by means of some specific FD 
approximations of the space derivatives y’ and y” [10]. Denoting the pulse 
magnitudes (in the space site 7) in the past, present and future time step by means 
of symbols p(i), pr(i) and f(i), there are obtained the basic KdV equation 
discretizations for the main studied y’ and y”” FD approximations (see below). 


2.1. The Landau-Paez discretization [11] 


Starting from the simplified (assuming V,, =0) expression of the KdV equation: 
OE SO a i al AO (2) 


the symmetric discretization: 


f() = POS [pri ++ pr@ + prG—Dh-LprG+)— prG—D)- 


— © [pr(i+2)+2pri—D) -2pr(i+ I) pr@-2)] (3) 
V5 


(where z and y are the FD time and space steps, respectively) is used. 
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2.2. The permanent 2 steps symmetric discretization [12] (pp. 271-275) 


The study of different FD discretizations pointed out that the asymmetric schemes 
lead quickly to the instability appearance. 


For this reason and simplicity, we studied previously [11] the permanent 2 steps 
symmetric FD discretization: 


f= p@)-y|pri+})- pr) |+a- pr(@)-[pré-1)=pr(i+))]+b-[pri-2)- pri+2)], (4) 


n-T aT ia 
with: OF a en OA DS (5) 
4 x x x 
where t and y have the same meaning (FD time and space steps, respectively). 
Given being these discretization schemes use the most simple (and rough) 
approximation of the first order derivative: 
j 1 : ‘ 1 ‘ : 
pf mek 2) = (+) zx, k-7]- yl) x, k-7}} = 7g Pret) pre NK) 


both these FD schemes lead to rather large simulation errors. 


2.3. The permanent 4 space steps symmetric discretization (this work) 
Starting from the substitution: 
u=V,+n-y, leading to: 4=n-y, w=n-y', u"=n-y", (7) 


one finds that the differential equation of the KdV pulses propagation is given by 
the simplified expression: 


u=—u-u'—d-u'''. (8) 


Additionally, instead of the 2 steps FD approximation (6) of the space first order 
derivative, it is used a 4 steps approximation: 


1 
yi yk Ose, ue = 2) -8y@=1) + 8yG +1) = yG + 2)}. (9) 


In this manner, it is obtained (using an averaged value of pr(i)) a FD schema using 
4 space steps: 


fi) = p@)- ia [pr +1) + pri) + pri-))]-[pri—2) -8pr(i—1) + 8pr(i +1) — pr +2)|- 
Xx 


d- 
— << [pri+2)+2pri-) -2prG +1) pri-2)] : a) 
Xx 
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Taking into account the expressions of the series expansions of the first order and 
of the third order derivatives by means of 2 space steps, and of 4 space steps, 
respectively (see e.g. [14], p. 162): 


1 - ' 
YG eet) = 5 G+ Dts kt MED x ke - AE BZ kt] ‘ (6’) 
and: YG ak t)= 1G 2) -8 yi NFB =yE+2)}4+ Sy (i-y,k-t), (7’) 
1 r V 
I to OQ HGS Di PG Dg VAIS NM 7), (1) 
y, 


as well as the rather small value of the space step y, it results that this FD schema 
(using a 4 steps symmetric discretization) is considerably more accurate (usually 
more than 10 times) than the previous 2 discretizations. 


3. Choice of the Uniqueness Parameters 

Besides the intrinsic uniqueness parameters of a Korteweg-de Vries pulse 
(soliton): the amplitude A (chosen as equal to 2), the 1-D speed inside an ideal 
(linear and conservative) medium V,, (chosen as 2), the non-linearity parameter n 
(chosen value 0.1) and the dissipation coefficient d (chosen value 0.2), we used 
the main following simulation uniqueness parameters: 

a) the FD space step y=dx (chosen value 0.15), b) the FD time step c=dr 


(measured in Vliegenthart [13] time units ty = [13], where — for the 


n-A+4d/ 7? 
above indicated values of the basic KdV parameters: ty =0.004195152), c) 
simulation space x,» x, d) initial position of the KdV soliton peak Xj, (usually 
chosen value in our simulations: x;,=200), e) number J of accomplished FD steps 


(iterations). An important simulation parameter - for the evaluation of the 
statistical parameters « skewness » and « kurtosis » of the simulated soliton: 


M M 
= oe (13) 
ml2? M2’ 
2 2 


where M,, is the statistical momentum of order 1 of the soliton components [10] - 
is the area A(/) under the soliton peak, equal initially to A;,; the change of this 
area after / steps (iterations) is denoted as: AA(/) [= A(/)—A;,,]. Another parameter 
whose values indicate the distortions magnitude in the /-th iteration, for the FD 
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site of index i, is the difference dy; (i- vy) between the simulated soliton component 
y,(i-y) and the true (physical) value 


yi“ 1-t)? NG N= EF D-YE- HTT). (13) 
Of course, for a global characterization of all distortions intervened after / 
simulation steps, the main parameter is the square mean error ¢ = (a7 G . ») of the 
simulated values relative to the true ones of the studied soliton components. 


The effective value of the KdV soliton speed is given by the expression : 


-A 
Vet. =Yoot—— 3 (14) 
for the above indicated numerical values, the effective speed is: Voge = 2.0666... 


4. Main obtained results and their interpretation 


4.1. The general time-space representation of the simulated soliton 


A general time-space representation of the simulated (by the numerical schema 
involving 4 space steps) soliton propagation is presented by Fig. 1. 


Fine K-VYoo=2; A=2: n=0.1; d=0.2: dt=1.6*tmx; dx=0.15 


(a.u.) 


1500 600 


erations \ 1000 


t/dt 0 0 Kidx 


400 


Fig. 1. Space-time representation of a KdV soliton propagation (the numerical values of the soliton 
shape and dynamic parameters, and those of the used FD simulation are written in this figure top). 
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In order to check both figure 1 and the above indicated theoretical model, we can 
estimate the soliton displacement along the OX axis during J = 2000 time steps: 
= Veff xIxdt — 2.066...x2000x1.6x 0.004195152 
dx 0.15 
agreement with figure 1. 


= 184.96 space steps, in 


4.2. Appearance of distortions 


Using the same values of the basic KdV parameters and of the FD simulation, we 
obtained the figure 2, involving the plot (in red) of the propagated soliton after the 
first 7 = 3) time step, and of the corresponding distortion 6y3(i-y). 


The amplitudes of distortions were multiplied by 10° to become visible. 


x 10° Fine K-Voo=2; A=2; n=0.1; 6=0.2; dt=1.6"tmx; dx=0.15 5% 10% Fine K-Voo=2; A=2; n=0.1, d-0.2; dt=1.6*tmx; dx=0.15 


(a.u.) 


k-u ; u/1000000 


A , 1 1 1 1 1 
in) 1 oo 0 300 700 500 600 tt} 100 200 300 400 500 600 
step j=d AS/S=-1.2532e-008% Wdx step 500 AS/S=-2.3075e-006% X/dx 


Fig. 2. The simulation begins with the step / =j =3, Fig. 3. The black plot indicates the distortions 


the first 2 iterations being injected analytically. (x 10°) for the iteration J = 500, while the 
We ask readers to observe that the distortions blue and red plots correspond to the initial 
begin with the first simulation iteration and present (J = 500) soliton positions, resp. 


4.3. Spreading and multiplication of distortions 


Figures 3 and 4 present the simulations of the KdV soliton propagation for the 
iterations J = 500 and J = 2000, respectively, as well as the corresponding 
distortions (multiplied by 10°), while figure 5 presents only the backward part of 
distortions for the iteration J = 2500. 


One finds that the pointed out distortions remain considerably weaker than the 
soliton components, but the number of distortions cycles (somewhat similar to a 
damped oscillation) increases considerably with the number of accomplished 
simulation iterations. 


Figures 6 and 7 present the same artifacts for both studied discretizations: the 
permanently 4 space steps discretization and the Landau-Paez one. 
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x10? Fine KVoo=2; A=2, n=0.1; d=0.2; dt=1.6%tmx; dx=0.15 x 10% Fine K-Voo=2; A=2; n=0.1; d=0.2; dt=1.6"tmx; dx=0.15 
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4 " 4 4 
100 150 200 300 
step =2500 
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step j=2000 
Fig. 4. For the iteration I = 2000, one finds both 
the soliton advance and the rise of distortions 
amplitude and number of cycles in the 


backward part of the soliton propagation. 


Fig. 5. For the iteration I = 2500, one finds that 
the distortions remain considerably smaller than 
the soliton amplitude, but the number of their 
(damped) cycles increase in the backward part. 


x10?  KVo0=2; A=2; n=0.1; d=0.2; dt=1.6%mx; dx=0.15 
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Fig. 6. Comparison of the distortions for the 
Landau-Paez (black plot) and of the permanent 4 
space steps discretization (red plot), respectively, 

corresponding to the iteration j = 100. 


Fig. 7. Comparison of the spreading and multi- 
plication of the distortions for the Landau-Paez 
and of the permanent 4 space steps discretiza- 
tion, respectively, for the iteration j = 2500. 


4.4. Dependencies of the relative square mean simulation error due to 
distortions on the number of accomplished iterations and the 
simulation space 


The above indicated results point out the increase of distortions amplitude with the 
number of accomplished FD iterations. Figures 8 and 9 present the plots of the relative 
square mean simulation error ¢ corresponding to different values of the time step 
t =dt , both for the 4 steps FD scheme and for the Landau-Paez scheme (see above). 
One finds: a) the considerably better accuracy of the 4 steps FD scheme relative to the 
Landau-Paez scheme (values of ¢ more than 10 times less for the 4 steps FD scheme), 
b) the abrupt rise (of burst type, leading to instability) of the relative square mean error 
é for values of the time step equal to (or larger) than 1.8 ty . 
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Fine K Voo=2; A=2; n=0.1; d=0.2; dx=0.15; Xmax=90 Landau K Voo=2; A=2; n=0.1; d=0.2; dx=0.15; Xmax=90 
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Fig. 8. Plots of the relative square mean error ¢ Fig. 9. The same plots in the same conditions, but 

versus the number of accomplished iterations for the Landau-Paez FD discretization schema. 

for different values of the time step T=dt in One finds that the values of the relative square 
Vliegenthart time f, and the 4 steps FD scheme. mean error are more than 10 times larger. 


Depending on the magnitude x», of the simulation space, the number of 


iterations corresponding to the launching of the « burst » phase of the distortions 
increases quickly (see Figures 10 and 11). In this phase, the square mean 
simulation error ¢ remains considerably less for the permanently 4 space steps 
discretization than for the Landau-Paez one. 


Fine Voo=2, A=2; n=0.1; d=0.2; dt=1.6"tmx; dx=0.15 


Landau K Voo=2; A=2; n=0.1; d=0.2; dt=1.6*tmx; dx=0.15 
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Fig. 10. Plots of the relative square mean error Fig. 11. The same plots in the same conditions, 
versus the number of accomplished iterations for for the Landau-Paez FD discretization schema. 
the permanently 4 space steps discretization and One finds that the values of the relative square 
different values of simulation space Imax =x,,,. mean error are of about 10 times larger. 


Finally, for very large values of the iterations number / = j (i.e. a huge «logical 
depth», which is also a Complexity feature), the plots of the relative square mean 
simulation error ¢ versus the iteration number reach a «saturation» phase. One 
finds (see figures 12, and 13) that in the frame of this «saturation» phase, the 
relative square mean simulation error tends to constant values of the magnitude 
order of 15+20%, for both studied discretizations. 
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As it is to be expected, the limit (« threshold ») values of the relative square mean 
simulation error decrease considerably (from approximately 35% to only 17.5%) 
for higher values of the simulation space (Imx = x,,, increasing values from 250 
to 1000 space steps), due to the reduction of the effects of interaction (reflection) 
with/on the boundary walls of the simulation space (see figures 12 and 13). 


Fine Voo=2; A=2; n=0.1; d=0.2; dt=1.6*tmx; dx=0.15 Landau K Voo=2; A=2; n=0.1; d=0.2; dt=1.6*tmx; dx=0.15 
40 . + T T + r r + r 


Imx=250 Irax=250 


imx= 1000 


it} 
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iteration j; (Nmx=10000) iteration j; (Nmx=10000) 


Fig. 12. Plots of the relative square mean error ¢ Fig. 13. The same plots in the same conditions, 

versus the number of accomplished iterations for for the Landau-Paez FD discretization schema. 
the permanently 4 space steps discretization and One finds that the values of the relative square 
different values of simulation space Imax =x,,,, mean error are almost the same as those for 


for the “saturation” phase (huge logical depth). the permanently 4 space steps discretization. 


One can find also that the numbers of iterations corresponding to the launch of the 
burst phase (figure 14) and the « threshold » (limit) values of the relative square 
mean simulation errors, reached in the « saturation » phase are practically the 
same for both studied FD discretizations. 


K Voo=2; A=2; n=0.1; d=0.2; dt=1.6%tmx; dx=0.15 K Voo=2; A=2; n=0.1; d=0.2; dt=1.6%tmx; dx=0.15 
r T r T T r 40 T T T r T r r T 
‘| | Fine ee 


Landau 


Imx=250 


0 Ooo 2000 3000 4000 5000 6000 7000 8000 9000 oo 1000 2000 3000 4000 S000 6000 7000 8000 9000 10000 
iteration j; (Nm»x=10000) iteration j; (Nm»x=10000) 


Fig. 14. Plots of the relative square mean errore Fig. 15. The same plots for both studied FD 
versus the number of accomplished iterations for discretization schemes, corresponding to the 
the permanently 4 space steps discretization and “saturation” phase. One finds that the limit 
for the Landau-Paez discretization scheme for the values of the relative mean errors are almost 

“burst” phase launch. the same for both these discretization schemes. 
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5. Study of the deformations of the simulated solitonic pulses shape 


In order to point out some specific deformations of the FD simulated solitonic 
pulses shape, there were calculated the moments of their components y,;(i) in the 


I-th iteration: 


imax F imax 3 
Mor= Lyr®, Mi= dLi-yi@, (15) 
imin imin 
and starting from the coordinate of « weight-center » of the simulated solitonic 
pulse components: (i), = /Mo » (16) 


the «centered» moments of the second, third and fourth order were also 
calculated: 


21 PS beedi) ,P-y1(0), Ma I (agile (0, Mar = SG PM, 17 


Pritt min tmin 


In this manner, it is possible to obtain the values of the skewness and kurtosis 
parameters for the /-th iteration: 


iy = A (12) 


While for the non-deformed (symmetric) launched KdV soliton: sg =0 and: 
kg =4.2, after several iterations both |s;| and |k; —ko| increase. As it is well- 


known, the non-null values of the skewness s correspond to asymmetric pulse, the 
reduced values of kurtosis k correspond to flattened pulses, etc. (see figures 16). 


¥y) 


Fig. 16. Comparison of the dotted plots corresponding to FD simulated KdV solitons with the 
continuous plots for the true (physical) KdV solitons, for different values of the skewness and 
kurtosis parameters, respectively. 
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K-Voo=2; A=2; n=0.1; d=0.2; dt=1.6*tmx; dx=0.15 


k; u; (au) 


426 426 430 432 434 436 438 
iteration j=2500 X/dx 
Fig. 17. Comparison of the plots for the Landau-Paez and permanently 4 space steps discretization 
schemes, with the red plot of true (physical) KdV soliton. 


Taking into account that for not too large numbers of iterations, the changes of 
skewness and kurtosis are rather small: |s;|<<1 and: |k; -ko|<«<1, figure 17 


presents only the top of the KdV simulated and true, resp. solitonic pulse in order 
to point out the changes of its shape and compare the results of the permanently 4 
space steps FD discretization with those of the Landau-Paez one. One finds again 
that the 4 space steps discretization leads to more accurate FD simulations. 


6. From the numerical artifacts to the numerical phenomena 


As it was shown already (see e.g. [4]), the detailed explanation of the 
mathematical mechanisms of the new found numerical artifacts (which transform 
these artifacts in numerical phenomena) is generally very intricate. As it can be 
expected, the corresponding difficulties are even higher for nonlinear physical 
systems. For this reason, we will examine here only the mechanism of the 
distortions appearance, beginning with the first iteration of the FD simulations. 


From the expression of the solution of the Korteweg-de Vries equation (1): 


y=A-cosh~? ® , where: @=4 x—[Voo tt }+| and: p=" 5 (15) 


it results that: y" =8A€° sinh of cosh ®—cosh™ 0) ; (16) 
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and: y” =-16AE sinh ®(2 cosh? @—30cosh™> ©+45cosh~/ ) . (17) 


The distortions amplitude 6y;(i-7)=y;(-y)-y@-yx,1-7) [defined by the expression 
(13)] involves the contributions of the corrections of the: (i) time derivative: 


r aes 
OV dy/ dt =e , (18) 


(ii) nonlinear term dy, , and of the : (ii1) dissipation term: 


2 


ova =. py" i. (19) 


While the expressions of the corrections of the time derivative and of the 
dissipation term, respectively, coincide both for the Landau-Paez and for the 
permanently 4 steps discretizations: 


ger73 
3 


3 
OY dy/ dt =— 4) -Asinh (3 cosh? © —cosh~/ o} (18’) 


(being of the magnitude order of 10°'° for the chosen input data, see paragraph 3), 
and: 


42 
byq = 12 4. Asinh @2cosh3 © -30cosh 9 +45 cosh~7 &) (19° 
mars (19°) 


(magnitude order 2x10~ for the chosen input data), the expressions, magnitude 
order and signs of the contribution dy, of the nonlinear term are different for the 
Landau-Paez and for the (permanently) 4 steps discretization scheme, 
respectively: 


B92 
OYLPn = — SEES. sinh ob cosh~> © —cosh~/ o) (20) 
(negative, of the magnitude order 4x10~° for the chosen input data), and: 
8ng ° yor neg 7 7 
Nan =— a Asinh ®(2 cosh ~ ®—30cosh ~ ®+45cosh " ®) , (20’) 


also negative but of the magnitude order of only 10'° 


One finds so that the dominant distortions term is the (negative) nonlinear term 
Oyppn (20) for the Landau-Paez discretization and the dissipation term dyg (19’) 
for the permanently 4 space steps discretization scheme. 
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From the above expressions, one finds that: 6dy,()=0 for sinh®=0, hence for: 
r=[YootS] 6 i.e. for the peak of the simulated soliton, in agreement with 


figures 2 - 7. Additionally, it results that at the right of the soliton’s peak, i.e. for 
®>0 we have: Sv4steps@) >0 and dyyp(i) <0 and — at the left of the soliton’s peak, 
i.e. for B<0: dy4steps@) <0, Svy_pp(i) >0, again in agreement with figures 2 - 7. 


Taking into account that for the Landau-Paez discretization the extreme values of 
distortions occur then for extreme values of y” , and the fourth order derivative of y is: 


yy =8Ae* (cost? &—15 cosh @ +15 cosh o), (21) 
it results that these extremes correspond approximately to the solutions of the 
equation: 15x? —15x+2=0, where: x= cosh 2 ®. 

One obtains: cosh® = esol ~1,09007, hence: e® ~1.09007+0.433892 and: 
154+ v105 


Dmx/ min = £0.4213165. Because for the used numerical values, the pseudo-wave- 
vector € (see equation (16c)) is: € ~ 0.2887, the relative (in respect with the KdV 
soliton’s peak) locations inyx; min Of the maximum and minimum of the Landau- 
Paez discretization distortion, respectively, are given by the relation: 


(a0) | 
: -[¥e i > , i — —maw min ~ 41 46005 , 
3 extr. 


[x—(Vo9 +nA/3)-t] 


hence Ainax/ min = A 


somewhat less than 10, again in agreement 


with Figure 6 and 7. 


As it concerns the position of the second “node” of the initial distortions (see 
figures 2 and 6), it is given for the Landau-Paez discretization by the condition: 


3cosh™> @-cosh=> ®=0, hence: cosh ®>,p =V/3, and: 
e?21P — 34/2 , obtaining: &y7p =1.146216. 


Finally, one obtains the relative (in respect with the KdV soliton peak location) 
position of the second “node” ( dy; (@ =0 ) of the initial distortion: 


Dorp 
&-dx 


Ai2nd nodesLP = (= 26.4684 for the used numerical values). 


One finds again the agreement with the simulation results presented by figure 6. 
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Similarly, for the permanent 4 space steps discretization, the position of the 
second “node” of the initial distortions is given by the condition: 


2cosh~* @—15cosh~* @+15cosh~° © = 0, 
leading to the equation: 2y? —30y+45 =0, where: y= cosh? ® ; (22) 
The solution of the equation (22) is: O2.=0.756676, leading to: 


® : 
Aiznd neath ZH (=17.4732 for the used numerical values). 
One finds that the first oscillation of the distortions of the permanent 4 space steps 
discretization is sensibly narrower than the corresponding one of the Landau-Paez 
discretization, again in agreement with figure 6. 


Of course, the accurate quantitative explanation of all pointed out numerical 
artifacts requires a rather difficult’ additional study. 


7. Conclusions 


A new permanently 4 FD space steps discretization was proposed and its 
numerical artifacts were systematically studied (for several uniqueness parameters 
of the studied solitons and simulations) and compared with those of the Landau- 
Paez «classical» discretization. Despite of the well-known difficulties of the 
numerical artifacts examination, the accomplished study succeeded to point out 
the possibility to explain quantitatively the mechanism of several numerical 
artifacts, which can be considered so as numerical phenomena. 


Taking into account the huge numbers of solitons applications [18], [19], and the 
efficiency of numerical simulations, the newly obtained results are useful for the 
improvement of the accuracy of FD simulations of the solitons propagation. 


'We will mention that many problems in the field of Numbers Theory are extremely difficult. E.g., 
the statement of the (Pierre de) Fermat's last (greatest) theorem was published in 1670 [15], by 
his eldest son — Clément Samuel Fermat, but its solution was found only in 1995 [16] by the 
American professor Andrew Wiles. Wiles describes ([17], p. 236) his experience of doing 
mathematics in terms of a journey through a dark unexplored mansion: “One enters the first room 
of the mansion and it’s dark. Completely dark. One stumbles around bumping into the furniture, 
but gradually you learn where each piece of furniture is. Finally, after six months or so, you find 
the light switch, you turn it on, and suddenly it’s all illuminated. You can see exactly where you 
were. Then you move in the next room and spend another six months in the dark. So each of these 
breakthroughs, while sometimes they’re momentary, sometimes over a period of one day or two, 
they are the culmination of, and couldn’t exist without, the many months of stumbling around in 
the dark that precede them”. 
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